/ mj t “\. ./ ^ .r *■ ^ .# #■ 


» 

^irnmr 

f 


/ ji V 


/^J 


; / 


\ 


OF 


KJ C I J. L 'w Li .i V •>» j C '.; \ V- .L V.. » ••’* '*• 

bOLLEGE 

ENGINEERING 


V n, 7 c 




VIRGINIA 
POLYTECHNIC 
INSTITUTE 

AND 

STATE 
UNIVERSITY, 


BLACKSBURG, 

VIRGINIA 




Final Report on NASA Grant No. NAG 3-593 

Thermodynamic Evaluation of Transonic Compressor Rotors 
Using the Finite Volume Approach 

for the period 
12/20/84 - 12/19/86 

by 

John Moore 

Professor of Mechanical Engineering 
Principal Investigator 

Stephen Nicholson 
Instructor 
and 

Joan G. Moore 
Research Associate 


Grantee Institution - 

NASA Lewis Research Center 
21000 Brookpark Road 
Cleveland, Ohio 44135 


Turbomachinery Research Group 
Report No. JM/87-4 

Mechanical Engineering Department 
Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 24061 



TABLE OF CONTENTS 

Page 

1. Abstract 3 


2. Development of an Explicit Time Marching Procedure 4 

for Laminar and Turbulent Flow 

- Summary Viewgraphs 

3. List of Project Reports and Papers 13 

Title Pages and Abstracts of Reports In Appendix A 

Copies of Papers In Appendices B and C 

4. Backflow - Extensions to the Computational Procedure 14 

Discretization of Convection Term 

Improved Pressure Interpolation for SBLI 
Evaluation of Turbulent Viscosity 

5. Backflow - Test Cases 22 

UTRC Separated and Reattached Turbulent Boundary Layer 

MDRL Diffuser - Strong Shock Case 

6. Mach Number Dependent Interpolation Formula 35 

for Density-Update Time-Marching Methods 

Appendices 

A. Title Pages and Abstracts of Reports 46 

B. An Explicit Finite-Voltme Time-Marching Procedure 59 

for Turbulent Flow Calculations 

« 

C. Explicit Finite-Volume Time-Marching Calculations 71 

of Total Temperature Distributions in Turbulent Flow 


- 2 - 


1. ABSTRACT 


In this final report, we summarize two years work developing 
computational capability to handle viscous flow with an explicit 
time— marching method based on the finite volume approach. For 
attached flow, our findings have been extensively documented, and 
our main object in this report is to present extensions to the 
computational procedure to allow the handling of shock induced 
separation and large regions of strong backflow. Two test cases 
are considered, the UTRC separated and reattached turbulent 
boundary layer and the strong shock case in the MDRL transonic 
diffuser G. The extended method has worked well on the UTRC flow 
with a boundary layer blockage of 58 percent and a maximum 
backflow velocity of 37 percent of the local maximum free-stream 
velocity. It has also worked well on the MDRL diffuser with a 
shock Mach number of 1,353 and a maximum backflow velocity of 
-71.7 m/s. 

A Mach number dependent interpolation formula for effective 
pressure has been developed for use in dens ity— update time- 
marching methods. This is a parallel development based on our 
earlier stability analysis which resulted in the M&M interpolation 
formula for effective density. 
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2. DEVELOPMENT OF AN EXPLICIT TIME MARCHING PROCEDURE 
FOR LAMINAR AND TURBULENT FLOW 

- SUMMARY VIEWGRAPHS 


- 4 - 


AN EXPLICIT FINITE-VOLUME TIME-MARCHING PROCEDURE 
FOR TURBULENT FLOW CALCULATIONS 

Start:: Denton explicit time-marching method. 
Allure - easy to understand method. 

Continuity 

dp + dpu + dpv + dpw = 0 
at dx dy dz 

6p = [- dpu - dpv - dpw ] 6t 
dx dy dz 

Momentum 

6(pu) = [steady eqn momentum error] 6t 
6(pv) = [ ] St 

6(pw) = [ ] St 
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Start: Denton explicit time-marching method. 

Questions : 

1. Is smoothing necessary for convergence of explicit method ? 


2. Why, at low Mach numbers, is the CFL criterion used to get the 
time step for the momentum equations ? 

( 6t = <Sx/ [velocity + speed of sound] ) 

3. Why not extend the method to laminar and turbulent flow ? 

What are the problems involved ? 

4. Why does he use an interpolated pressure in the momentum 
equation for transonic flow ? 

Can we show why and when it is stable ? 

5. How can the method bo extended to separated flow ? 



ANSWER > 

Development of Explicit method for calculation of 
Inviscid, Laminar or Turbulent Flow 
Mach number = 0 to >2.5, including shocks 
Economical - grid points 
With or without separation 
Tested on 2-d duct flows 
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1. Is smoothing necessary for convergence of explicit method ? 
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( "New" control volume, traditionally used for boundary layers ) 
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2. Why, at low Mach numbers, is the CFL criterion used 
to get the time step for the momentum equations ? 

( 6t = 6x/ [velocity + speed of sound] ) 

CONSERVATIVE FORM OF MOMENTUM EQUATION 

d(pu) + V'pji u = - dp + . . . 
at ax 

U [ ap + V*pji ] + pdu + PlA’Vu = -dp + . . . 

at at ax 

continuity 

included, therefore 6t = <St _ 

coQu mom 

stability analysis, continuity and momentum > 

CFL condition 6t = 6x/(u+c) 

CONVECTIVE FORM OF EQUATION 

pau + pu’Vu = -dp + . . . 
at dx 

Stability analysis, momentum equation > 6t = 6x/u 
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3. Why not extend the method to laminar and turbulent flow ? 
What are the problems involved ? 

RESOLUTION OF TRANSVERSE PRESSURE GRADIENT 



Flat plate turbulent boundary layer dp/ay n 0 

P = pRT 

6p = RT 6p - pM 6U 

20p 

6p dependent on continuity and momentum errors - 
stability is highly grid and 6t dependent 
difficult without smoothing 

Borrow idea from pressure correction methods - 
6p depends ohly on continuity error. 

6p = RT 6p 

= [ -dpu - apv - dpw ] RT 6t 
6x ay az 

Stable without smoothing. Multi-volume approach 
needed for highly nonuniform 6y grid spacing. 
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4. Why does Denton use an interpolated pressure in 
the momentum equation for transonic flow ? 

Csuii we show why and when it is stable ? 


Wamt apu + apv + apw = 0 

ax ay az 

1-D stability analysis. p = p + 6p 

u = u + 6u 
6u r. ~ C a<6p)/dx 
<5p ^ RT ap 

Interpolated pressure 

p^=(pRT)^ or Pj^=.5[(pRT)^^j^+(pRT)^_j^] or . . . 


> 


Ai«P. 


= continuity error, 


+ . . . 

each control volume 


Explicit method approximation 

[l/6t] = continuity error for control volume 

Stability requires positive and dominant. 
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5. How can the method be extended to separated flow ? 
(a) UPWIND DIFFERENCING 


> u 

O 

• 

o 


li-l 

J 

1 i 

i+1 

u > 0 

U du/dx = 



u 


u < 0 

U du/dx = 

(-U) (u^ 



> positive coefficient for 

(b) UPWINDED CONTROL VOLUMES 

- control volumes depend on local u 


> positive coefficient for 
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4. BACKFLOW - EXTENSIONS TO THE COMPUTATIONAL PROCEDURE 
4a. Discretization of Convection Terms 


The momentum and energy equations are discretized over control 
volumes fixed relative to the grid points. Central differencing is 
used except in regions where there are large cross flows or 
backflow. In these regions a side upwind or reverse upwind 
differencing is used for stability. The details are as follows. 

Control volume for momentum or energy for point i+l,j. 

X X j + i 

N 


W X X E j 


S 

X X j-1 

i i+1 


Bulk flow 
direction 
> 


Convection of property where 

(^ = u for x-momentum 

(^ = V for y-momentum 

<i> = h for energy (enthalpy) equation. 

Convection term integrated over control volume 

dVol = /V'puit dVol - ^j^/V’pu dVol 

= ^<l)Pli*dA - 4)j(j/Pii*dA 

= + (pU’A> 3<^3 + (Pii*A)g<l>g + (Pil*A)(^4>j^ 

- <l>j4[(Pil*A)jj + (pu*A) 3 + (PU*A)g + (Pli*A)j(] 

We wish to express this in terms of the (^*s at the grid points, 
i.e we wauit the equation in the form 


/pu*7ij) dVol 


Vi+1, j Vi, j ‘^EE*i+2,j 


The coefficients C are determined from the mass fluxes through the 
sides (pu*A) and the discretization choice for (Jig, (i>g, and 


For stability we wish the center point coefficient, C„, to be 

hi 

positive auid greater than the sum of the other positive 
coefficients . 



When [(Pii'A)jj+(PU’A) 5 +(PU'A)g+{Pii*A)^] > 0 
take (J)j^ = ^ 

to give a negative contribution to 
When [(Pii-A)j,+ (PU-A)s+(PU-A)E+(Pii-A)pj] < 0 
take = *i+i, j 

to give a positive contribution to Cg. 

<J»g and <t)^ 

When (pu*A)g > 0 
= ‘**i+l,j 

This centered evaluation of (second order accurate) gives a 
positive contribution to Cg. 

When (Pii’A)g < 0 
" ^i+2,j 

This upwind evaluation of (|ig (first order accurate) gives a 
negative contribution to Cgg. 

This also determines (|>^ since (tg for one control volume is (t>j^ for 
the next control volume. 
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and <i>g 


The geometrical centered evaluation for is 

where F is the fraction of the distance of the North face between 
the grid points. 


j+1 


-N - 


I b 


F = a/b 


For accuracy this centered evaluation should be used whenever 
possible. 

For stability when 


take 


> 0 and (pu*A)g > 0 


F ^ (Pii*A)E/(PiA‘A)jj 


When the inequality is chosen, which for equal grid spacing will 
occur when (pu’A)jj > 2(pii*A)g. 

(pu-A)^«l>jj = j + j) t(Pil’A)jj - (PU’A)g] 

When (pu*A)jj > 0 and (Pii’A)g ^ 0 

(the primary flow is locally backwards or zero relative to the 
bulk flow direction) 


‘*’N " “^i+l, j 


take 

By symmetry for {pu'^)g > 0 and (PU'A)g > 0 
take F 2 1 - ( PM’A)g/(PlA'A )5 

For (pu’A)g > 0 and (Pii*A)g S. 0 




take 




When (pu*A)jj is negative, (pii’A)g is positive for the next, the 
j + 1, control volume, so that is determined from <|)g for the j+1 
control volume. Similarly when (pu"A)g is negative, <|)g is 
determined from for the j-1 control volume. 

Comparison with earlier scheme. 

In these terms, Nicholson (Section 3, Report 5, JM/86-6) 
considered only positive values for (pu'A)g, i-e. no reverse flow. 

The formulae he used for and <1)^^ were the seune as given 

here. However the upwinding he took for the cross flows was 
different. In particular when (pu*A)jj was positive, <|)jj was 

evaluated using 

F i 1 i.l)/(PU-4)N + 1 

Taking F A ( Pii'A)g/( gives lower and hence more conservative 
values of F when (pii‘A)j^ > 2{Pii‘A)gp i.e., where for uniform grid 
spacing, the geometric F may not be used. 



high low 

cross flow cross flow 


F used in calculations. 
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Table 1 SIDE UPWINDING ALGORITHM 
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4b. Improved Pressure Interpolation for SBLI 

The same Mach number dependent pressure interpolation formula 
for the calculation of the density used earlier is also used here. 
However to converge calculations with strong shock boundary layer 
interactions, i.e. with shock induced separations, the Mach number 
used in the formula needed to be changed from the local Mach 
number to the local free stream Mach number. Since stability is 
compromised if the Mach number is underestimated but not if it is 
overestimated, for simplicity, the value used was the largest Mach 
number on the relevant pair of i-surfaces. 


X 

1 

X 

X 

1 

X 

X 

X 

1 

1 

X 

X 

1 

1 

X 

X 

X 

1 

X 

X 

1 

X 

X 

X 

1 

X 

X D - 

1 

X 

X 





1 



X 

1 

X 

X 

1 

X 

X 


M for these points 
max 


"iti = 5^ ‘Pi ^ <Pi+i-Pi-i> 

Sjj(M) + »j(M) + SjfM) = 1 
M = Mach number = maximum Mach number at 


" <Pi+l-Pi-2>l 

planes i and i+1. 
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4c. Evaluation of Turbulent Viscosity for the Present Test Cases 


The turbulence model used in the calculations is a Prandtl 
mixing length formulation 

= pL^"du/dy" 

where the mixing length L is the smaller of 

0. 41y (with a Van Driest correction) 
or 0.086. 

Sometimes it is difficult to determine an appropriate boundary 
layer thickness, 6. 

For the present test cases 

(a) UTRC, boundary layer separation geometrically triggered and 

(b) Sajben diffuser, separation induced by shock, 

the boundary layer thickness used to calculate L was ch 2 uiged to 
obtain a reasonable separation when compared with the 
measurements. The details of what was used to determine the 
effective boundary layer thickness follow. 

Determining the edge of the boundary layer. 

The location of the edge of the boundary layer is determined 
by the total pressure gradient |dp^/dyl. In particular a search 

starts from outside the boundary layer (in the middle of the duct) 
and proceeds towards the wall to locate where 

dp, _ (p_ - P) *DPFACT 

Li ~ — u ! 

dy (local duct width) 

This is the edge of the boundary layer for the mixing length 
calculation. The larger the DPFACT the smaller the boundary layer 
thickness. 


Case DPFACT 

(a) UTRC, sep. b. 1. zTo 

(b) Sajben, P_/P_ = 0.722 2.5 

6 O 
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"Time” lag for boundary layer thickness. 


For case (b), Sajben P^/Pq = 0.722, shrinking 6 by increasing 

DPFACT was insufficient to correctly obtain the separation induced 
by the shock. Qualitatively since turbulence is convected with the 
flow there needs to be time for the turbulence to change - it does 
not change suddenly. This was qualitatively introduced into the 
calculation by lagging the boundary layer thickness used for the 
calculation of L by 5 grid points. (The lag is between 0.5 and 
0.75 throat heights through the separation region. ) In particular 
after i = 40 (x/h=1.7, upstream of the shock at x/h~2.4 but well 
downstream of the throat at x/h=0) the mixing length in the outer 
part of the boundary layer was obtained using 

L(i) = 0.086(i-5). 

The time lag was used only for case (b). 
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5. BACKFLOW - TEST CASES 


The extensions to the computational procedure described in 
section 4 were necessary for modelling two extreme cases of 
separated flow, the UTRC separated and reattached flat plate 
turbulent boundary layer (NASA Contract NAS3-22770, reference 1) 
and the MDRL transonic diffuser flow with a strong shock (MDRL 
Report No, 81-07, reference 2). These cases exhibit large boundary 
layer blockage (displacement thickness/local duct height), large 
backflow velocities, relative to the free stream velocity, and 
high rms/meaui turbulence levels in the backflow region. The 
maximum boundary layer blockages were 58 percent (fifty eight!) in 
the UTRC low speed = 27 m/s) flow and 27 percent in the MDRL 

diffuser with a shock Mach number of 1.353. The maximum backflow 
velocities were 37 percent and 25 percent, respectively, of the 
local maximum free-stream velocity. The ratios of rms/mean axial 
velocities at the locations of maximum reverse flow velocity were 
35 percent in the UTRC flow and 66 percent in the MDRL diffuser. 

If the backflows in the two cases were varying sinusoidally, these 
values would correspond to maximum backflow velocities of -5.4 + 
2,7 m/s amd -71.7 + 67.3 m/s, respectively. 

UTRC Separated and Reattached Turbulent Boundary Layer 


The geometry and streamlines for flow through the UTRC test 
section are shown in Fig. 1; and laser doppler velocity 
measurements are shown in Fig. 2. The corresponding calculated 
velocity vectors together with the locus of points for which U=0 
are seen in Fig. 3. The size of the reverse flow region is well 
modelled, and the maximum calculated reverse flow velocity of -4.1 
m/s agrees well with the measured maximum value of -5.4 m/s. This 
good agreement for the reverse flow leads to reasonable agreement 
between the measured and calculated values of skin friction 
coefficient in the separation zone, as shown in Fig. 4. The 
calculated locations of separation and reattachment are seen to be 
close to the measured locations. The good agreement in the 
separated flow region was obtained with the Premdtl mixing length 
model by reducing the turbulent viscosity in the boundary layer as 
discussed in section 4 (i.e. by using DPFACT = 2.0). This then 
gave a corresponding decrease in the calculated skin friction 
upstream and downstream of the separation zone, as seen in Fig. 4. 
We conclude that the present explicit ccmputational procedure can 
be used for flows with extensive and strong backflow but that a 
more sophisticated turbulence model is required. 
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MDRL Diffuser - Strong Shock Case 


With a back pressure, Pg^it^^o inlet' 0.722, the MDRL diffuser 

G had a shock Mach number of 1.353. Shock induced separation 
occurred in the turbulent boundary layer on the curved top wall. 
This contrasts with the case of 0.805 pressure ratio which gave a 
shock Mach number of 1.235 and no separation. In this section, 
results of calculations for these two flows will be compared, with 
particular attention being given to the backflow in the strong 
shock case. 

The calculated shock locations are clearly seen for the two 
cases in the contours of static pressure in Fig. 5. The strong 
shock is located further downstream and shows evidence of a lambda 
foot at the curved top wall. The computed shocks are both quite 
sharp as a result of the use of the MStM pressure interpolation 
formula (see section 3 of this report, reference 3). 

For the strong shock case, the computed and measured static 
pressure distributions on the top wall are compared in Fig. 6. The 
computed shock is just downstream of the measured location and is 
therefore somewhat stronger with a shock Mach number of 1.39. 
Upstream of the shock the static pressures are indistinguishable; 
but downstream the calculated static pressures are consistently 
higher than those measured, perhaps partly because of 
three-dimensionality in the measured flow. 

The Mach number contours in Fig. 7 show the flow accelerating 
up to the shock and decelerating downstream. The top wall boundary 
layer thickens appreciably more through the strong shock. This is 
seen also in the velocity vectors of Fig. 8, which show the 
separation bubble downstream of the strong shock. Fig. 9 shows 
this calculated backflow in more detail, and for comparison the 
magnitude and possible variations of the measured backflow are 
also shown. The maximum calculated backflow velocity of -87.7 m/s 
agrees quite well with the maximum measured value of -71.7 m/s. 

Figures 8 and 9 demonstrate quite graphically the significant 
blockage caused by the separation bubble; and this is also seen in 
Fig. 10, which shows contours of total pressure. 

We conclude that calculations of diffuser flows with strong 
shocks and shock induced separation can be performed with the 
present explicit method. As discussed in section 4, this 
calculation required a time lag of the turbulent viscosity to give 
a reduced viscosity in the separation bubble. In fact, this simple 
modification to the turbulence model produced a dramatic upstream 
movement of the shock and the calculation rapidly converged on a 
shock location close to that measured. Again this suggests the 
need for a more sophisticated turbulence model. But the present 
study of strong backflows has clearly demonstrated that they can 
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be modelled with an explicit method based on the finite volume 
approach . 
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Calculated velocity vectors and the locus of points for which U 
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Skin Friction Coefficient on Test Surface (z=0 , after Patrick [1]). 
Comparison of calculation results and measurements. 
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Fig. 6 Computed and measured static pressure distributions 
on the curved top wall of the MDRL diffuser. 

Strong shock case; Pgxit'^^Oinlet “ ^.722. 


- 30 - 




Fig. 7 SAJBEN TRANSONIC DIFFUSER Mach Number, Contour Interval = 0.05 



0.805 




Calculated velocity vectors upstream of and within the separation bubble 





Fig. 10 SAJBEN TRANSONIC DIFFUSER Total Pressure > Pr/Pn^ Contour Interval = 0.025 


6. MACH NUMBER DEPENDENT INTERPOLATION FORMULA 
FOR DENSITY-UPDATE TIME-MARCHING METHODS 


A 1-d stability analysis of density-pressure relations used in 
the computation of transonic flow was performed in Report No. 
JM/85-11 (see section 3 of this report, reference 3). Here we give 
a parallel development of a density interpolation equation for 
effective pressure for use in density-update methods. The formulae 
considered are tested using the density-update scheme outlined in 
Table 1. 

Downwind Effective Pressure 


In section 2.8 of reference 3, we considered an inconsistency 
in the pressure-density relation such that the pressure used in 
the momentum equation is offset by one grid point from the density 
used in the continuity equation, i.e. 

Pi = Pi+iRT. (40) 

In a density update method this may be viewed as an effective 
pressure evaluated downwind of its point of use in the momentum 
equations. This pressure-density relation was found to be stable 
for all Mach numbers, but it results in poor shock capturing as 
the calculated shock is spread over numerous grid points. Fig. 1 
shows the calculated and theoretical pressure distributions for a 
1-d calculation with a nominal shock number of 1.45. 

Mach Number Dependent Interpolation Formula for Effective Pressure 


In section : 
number is high, 
gas equation of 

P 


.6 of reference 3, 
the density update 
state satisfied at 


we saw that when 
method is stable 
each grid point. 


the Mach 

with the ideal 

i.e., 

(73) 


Since this is the correct pressure-density relation for ideal 
gases it should be used where feasible. In this section, we will 
start with a generalized density interpolation equation for 
effective pressure 






(74) 


with 


e.RT. 

= P 1 1 


(75) 
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and 


1 


(76) 


for second order accuracy. 

We seek Mach number limitations to a^, a^ and ag using the 
stability criterion that 

the center point coefficient must be greater than 
the sum of the other positive coefficients, 

^°®^center ^ 


Substituting 



2 3 




- Vi +3 


(77) 


into Eq. 25 of reference 3 and rearranging in terms of the 
coefficients of each 6p^, a^, a^^, and ag, yields (neglecting 

variations of A, u and c with i) 


Ac 

«(M+1> 


{ ( 
( 
( 


i ^ ^^i+4 


1 ai - I a2 ) «Pi^3 


-1 


+ a_ - a. + 1 a« ) 6p 


i+2 


( 2+M«(M+l) - 3a 


- i »2 > ■*Pi+l 


(-1-MS(M+1) + 3a + a, + 2 a, ) «P, 

3 


( 


- - 1 a, - 1 a« ) 6p. 


2 S “2 


i-1 


= A 


change, i ' 


(78) 


Now let us consider a simple second order scheme with a 2 = 0 and 
a^ = 1 - a^, and find limiting values of a^. From Eq. 74, it is 
obvious that we should consider only values in the range 
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( 79 ) 


0 i. ^ 1. 
o 


In this range, the coefficient of is positive or zero, and 

for the coefficient of (the center point) to be greater than 


the coeffcient of we require 


a^ < S +2 M«(M+1) . 

° 5 5 

The coefficient of 6p^ is positive when 

2a - M«(M+1) > 0 or M(M+1) < 2a /«. 

o o 

In this region, from Eq. 46, we require 

2 + M«(M+1) -3a >l-a +2a - M«(M+1) 

o 9 _Q o 

2 

or a^ < 1 + 4 MS(M+1) . 

° 3 9 


(80) 


(81) 


(82) 

(83) 


This criterion is more restrictive that Eq. 80 and the 
corresponding stability limit is shown as a function of Mach 
number in Fig. 2 for the conservative case of < = 1.0. 

A set of equations for a^, a^^ and agj which satisfy the 

stability criteria (Eqs. 80 and 83) and give second order accurate 
interpolation (Eq. 76) has been selected; that is 

a = 4 M(M+1) 

° 9 

a. = 1 - a (84) 

1 o 


This Mach number dependent formulation for a^ and a^^ is shown in 

Fig. 3. These equations are referred to as the M&M Mach number 
dependent interpolation formula for density-update time-marching 
methods . 
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Computational Tests of the M&M Density Interpolation Formula 


In this section, results of shock capturing with the MStM 
formula (Eq. 84) in the density update method (Table 1) are 
presented for Denton’s 1-d nozzle. 

Calcula tion Details 

Number of axial grid points 
At inlet i 

For air k 

^exif^^t, inlet 
Results 

The variations of static pressure, Mach number, and total 
pressure for all three back pressures are shown in Fig. 4. All 
three shocks are captured over four steps. The upstream side of 
the shock is sharply defined with only minor deviations from the 
theoretical 1-d solution. On the downstream side, there is a small 
overshoot and undershoot in static pressure and Mach number ever 
two steps; the total pressure distributions show no overshoots or 
undershoots and show a sharp decrease over two steps. 

Concluding Remarks 

It is hoped that the M&M density interpolation formula will be 
useful to those organizations like NASA Lewis who are using 
density-update time-marching codes. It is also hoped that the 
stability analysis performed under this NASA Grant will be 
enlightening to users and developers of time-marching codes. 


= 46, 6x = 1 

= 1, M = 0.80 

=1.4, R = 287. J/kgK 

= 0.85, 0.80, 0.75 
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Table 1. 


Outline of Effective Pressure Method with Different Time Steps 
Density Update Scheme. 

UNKNOWNS ( 2 - DIMENSIONS) 

^tU,v,y?u) ,y»v) , ^e^),h^,P,T 

CONTINUITY 

1^ + - 0 Cl) 

%f> - (-<^y;u)-dA) ^tyVol 

f 

P - ^RT 

MOMENTUM (INVISCID) 
h (^u) + y»u) u = -V p 


USING EQ. 2 - (u or v) TIMES EQ. 1 

!>U - (-§ (uy?u)«d4+ pi^»dA) 

+ (j}\^»dA))*^ t^/Vol 

■ (-^(v^u)«dA + pjL«dA) 

+ y>u)«dA))*S 


u ■ u + S u 

V ■ V +Sv 


(^) - y?u 

(y>v) -y>v 


ENERGI 


h 

o 

T 


* constant 

" (h^ - u^+v^) /C 
0 — = — P 
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STBBILITT LIMITS. RO+Rl'l 



0.00 0.20 O.UO O.SO O.BO t.OO 1.20 I.MO l.BO 1.00 2.00 2.20 2.00 2.60 2.00 3.00 

MflCH 

Fig. 2 M & M density interpolation formula for use with 

Effective Pressure, Density Update, Time Marching Methods. 

Stability limits for a^ when a 2 - 0 and a^ + a^^ - 1. 
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MfiCH NUMBER DEPENDENT fl*S WITH RO+RUl 



0.00 0.20 0.<0 0.80 0.80 1.00 1.20 I.UO 1.80 1.00 2.00 2.20 2.80 2.80 2.00 3.00 

MRCH 


Fig. 3 M & M Mach number dependent values for the density interpolation 

coefficients, a^ and in the Effective Pressure Method. 
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DENTON ID EXAMPLE ftO=Mx (M+1) x4/9 = 


pw = 


p/p 


t, Inlet 



o 


0.0 10.0 20.0 30.0 > 40.0 50.0 


X 


^^8* Calculated 1— D solution for Denton's nozzle 

using M & M density interpolation formula 
with the Nlcholson/Moore Effective Pressure, 

Density Update Method. 

Calculations for three exit static pressures at x = 46., 
inlet ' “•*5. 0.80, and 0.75. 
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DENTON ID EXAMPLE 


RO=Mm CH+1) mI4/9 R1 = 1-R0 



X 


Fig. 4c PTOT 
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DEVELOPMENT OF A FINITE VOLUME TIME MARCHING METHOD 


ABSTRACT 

The Objective of the current work is to develop and demonstrate a 
Navier— Stokes approach for transonic flow which includes viscous 
terms in the finite-volume method. The accuracy of the 
computational method will be verified using a transonic diffuser 
as a test case. The computational goal is to calculate the flow 
in sufficient detail and with sufficient accuracy that the loss 
generating mechanisms can be studied to assess the sources of 
inefficiency in the transonic diffuser. The purpose of this 
report is to document progress made in the development of the 
time-marching finite-volume method from September 1984 to December 
1984. 
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EXTENSION OF A FINITE VOLUME EXPLICIT TIME MARCHING METHOD TO 
LAMINAR AND TURBULENT FLOW 


ABSTRACT 


This report documents progress made in extending the finite 
volume explicit time marching method to laminar and turbulent 
flow during the time period from January to May 1985, The work 
done is under NASA grant NAG 3-593. Previously, extensions had 
been made to the finite volume method to improve the accuracy of 
the calculation of total pressure in compressible invlscid flow. 
These changes are documented in reference 1 . The current work 
extends these ideas and develops new ideas which allow the 
calculation of laminar and turbulent boundary layers in internal 
flows. The method is verified using four test cases with 
free-stream Mach numbers ranging from .075 to 1.20. 
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Annual Report on NASA Grant No. NAG 3-593 

Thermodynamic Evaluation of Transonic Compressor Rotors 
Using the Finite Volume Approach 
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12/20/84 - 12/19/85 

by 

John Moore 

Professor of Mechanical Engineering 
Principal Investigator 

Stephen Nicholson 
Instructor 
and 

Joan G. Moore 
Research Associate 


Grantee Institution - 

NASA Lewis Research Center 
21000 Brookpark Road 
Cleveland, Ohio 44135 


Turbomachinery Research Group 
Report No. JM/85-11 

Mechanical Engineering Department 
Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 24061 


Abs tract 


Summer research at NASA Lewis Research Center gave the opportunity 
to Incorporate new control volumes In the Denton 3-D finite-volume tlme- 
marchlng code. For duct flows, the new control volumes require no 
transverse smoothing and this allows calculations with large transverse 
gradients In properties without significant numerical total pressure 
losses. 

The summer research also pointed to possibilities for Improving the 
Denton code to obtain better distributions of properties through 
shocks. Much better total pressure distributions through shocks are 
obtained when the Interpolated effective pressure, needed to stabilize 
the solution procedure. Is used to calculate the total pressure. This 
simple change largely eliminates the undershoot In total pressure down- 
stream of a shock. Overshoots and undershoots In total pressure can 
then be further reduced by a factor of 10 by adopting the effective 
density method, developed at VPI&SU, rather than the effective pressure 
method. Use of a Mach number dependent Interpolation scheme for pres- 
sure then removes the overshoot In static pressure downstream of a 
shock. 


The stability of Interpolation schemes used for the calculation of 
effective density Is analyzed and a Mach number dependent scheme, the 
M&M formula. Is developed. This formula combines the advantages of the 
correct perfect gas equation for subsonic flow with the stability of 2- 
polnt and 3-polnt Interpolation schemes for supersonic flow. 
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ABSTRACT 


This report documents progress made in refining and improving the finite*volume explicit time 
marching method (1, 2, and 3 ) during the time period fi-om January to May 1986. The work is 
done under NASA grant NAG 3*593. Previously, extension had been made to the finite volume 
method to 

1. improve the accuracy of the calculation of total pressure in inviscid flow (1). 

2. extend the method to allow the calculation of laminar and turbulent boundary layers in internal 
flows (2). 

3. improve the shock capturing properties of the method by introducing a Mach number de* 
pendent interpolation scheme for the pressure used in the calculating the density (3). 

The current work extends these developments by 

1. using the new pressure interpolation scheme in two dimensional viscous calculations. 

2. including a more complete description of the viscous stresses. 

3. introducing a criteria for the transverse upwind differencing which is a fiinction of the ratio of 
transverse and streamwise mass fluxes. 

4. allowing the calculation of internal flow where boundary layers are present on both wall of the 
duct. 

Specifically, this report is broken up into three sections. Section 1 discusses in detail the maimer 
in which the viscous stresses are evaluated in the non*orthogonal, non*uniform grid. Section 2 in- 
vestigates the convergence and presents results for calculations of laming flow in a converging duct. 


- 54 - 



Section 3 presents results for calculations of transonic turbulent flow in a converging<diverging 
nozzle; the results are compared with Sajben's measurements and calculations by other authors. 
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Extension of the Finite Volume Method 
to Laminar and Turbulent Flow 
by 

Stephen Nicholson 
John Mooie, Chairman 
Mechanical Engineering 
(ABSTRACT) 

A method has been developed which calculates two-dimensional, transonic, viscous flow in ducts. 
The finite volume, time marching formulation is used to obtain steady flow solutions of the 
Reynolds-averaged form of the Navier Stokes equations. The entire calculation is performed in the 
physical domain. The method is currently limited to the calculation of attached flows. 

The features of the current method can be summarized as follows. Control volumes are chosen so 
that smoothing of flow properties, typically required for stability, is not needed. Different time steps 
are used in the different governing equations to improve the convergence speed of the viscous cal- 
culations. A new pressure interpolation scheme is introduced which improves the shock capturing 
ability of the method. A multi-volume method for pressure changes in the botmdary layer allows 
calculations which use very long and thin control volumes ( length/height S 1000). A special 

discretization technique is also used to stabilize these calculations which use long and thin control 

\ 

volumes. A special formulation of the energy equation is used to provide improved transient be- 
havior of solutions which use the full energy equation. 

The method is then compared with a wide variety of test cases. The fieestream Mach numbers 
range from 0.07S to 2.8 in the calculations. Transonic viscous flow in a converging diverging nozzle 
is calc ulat ed with the method; the Mach number upstream of the shock is approximately 1.25. The 
agreement between the calculated and measured shock strength and total pressure losses is good. 
Essentially incompressible turbulent boundary layer flow in an adverse presstire gradient is calcu- 
lated and the computed distribution of mean velocity and shear stress are in good agreement with 
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the measurements. At the other end of the Mach number range, a flat plate turbulent boundary 
layer with a ireestream Mach number of 2.8 is calculated using the full energy equation; the com- 
puted total temperature distribution and recovery factor agree well with the measurements when a 
variable Prandtl number is used through the boundary layer. 
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An Explicit Flnlte-Voluae Tlac'-Harchlng Procedure 
for Turbulent Flow Calculations 

Stephen Nicholson, Joan G« Moore and John Moore 

Mechanical Engineering Department 

Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 24061 

1 . SUMMARY 

A method has been developed which calculates two- 
dimensional, transonic, viscous flow In ducts. The finite- 
volume, tlme-marchlng formulation Is used to obtain steady 
flow solutions of the Reynolds-averaged form of the Navler 
Stokes equations. The entire calculation Is performed In the 
physical domain. 

The features of the current method can be summarized as 
follows. Control volumes are chosen so that smoothing of 
flow properties, typically required for stability. Is not 
needed. Different time steps are used In the different 
governing equations. A new pressure Interpolation scheme Is 
Introduced which Improves the shock capturing ability of the 
method. A multi-volume method for pressure changes In the 
boundary layer allows calculations which use very long’ and 
thin control volumes (length/height - 1000). The method Is 
then compared here with two test cases. Essentially Incom- 
pressible turbulent boundary layer flow In an adverse pres- 
sure gradient Is calculated and the computed distributions of 
mean velocity and shear stress are In good agreement with the 
measurements. Transonic viscous flow In a converging diver- 
ging nozzle Is calculated; the Mach number upstream of the 
shock Is approximately 1.25. The agreement between the 
calculated and measured shock strength and total pressure 
losses Is good. 

2. INTRODUCTION 

The finite volume method has been used extensively to 
solve the Euler equations for transonic flow Including flow 
at high Mach numbers. In Internal aerodynamics, McDonald [1] 
was the first Investigator to use the time marching finite 
volume method. Denton [2] extended McDonald's finite-volume 
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method to three dimensions. Versions of Denton's method have 
been used In lnvlscld~vlscous Interaction programs for 
turbomachinery calculations [3-5], 

The scope of the present work was to extend a finite 
volume method like that of Denton's to be able to calculate 
laminar or turbulent flow In ducts. The new method has the 

capability to calculate subsonic as well as transonic flow. 

3. GOVERNING EQUATIONS 

The unsteady foirra of the continuity equation, the x- 
momentum equation, and the y-momentum equation. In Integral 
form, are used to obtain a steady-state solution for flow 
through 2-dlmenslonal ducts. The Ideal gas equation of 
state, the assumption of constant total temperature, and a 
Prandtl mixing length turbulence model complete the governing 
equations needed to solve for the unknown variables p, u, v, 
P, y, and T. 

For a finite control volume where we can assign one 
value of density to the control volume, and for a finite time 
step, 6t, continuity states that, 

- p" = (5p » -[// pu • dA] (1) 


where the Integral Is evaluated explicitly at the current 
time step, n. In arriving at an expression which relates the 
pressure change directly to the continuity error, we will 
assume that changes In temperature are small In comparison to 
other changes for one time step. Thus, we can relate changes 
In pressure to changes In density through the Ideal gas 
equation of state, 

pn+1 . pn . 5 p . ^6^ ^ 2 ) 

For the method Introduced In the current work, a non-conser- 
vative form of the unsteady momentum equation Is used. The 
non-conservative form Is used because It allows the use of 
different time steps for the continuity and momentum equa- 
tions. The differences between the non-conservative and 

conservative forms of the unsteady momentum equations are 
associated with the unsteady and convective terras. Speci- 
fically, we note that 

3(pu) 3u 

+ V • pu u » p .jZ + pu . Vu (3) 
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and the right hand side of Eq. (3) can be rewritten as 


3u 

p ^ + pu • Vu 


3u 


pu u - u(V 


pu) 


( 4 ) 


When the right hand side of Eq. (4) is combined with the 
pressure and viscous terms, the momentum equation in integral 
form becomes 


(u)”"*"^ - (u)’^ * 4(u) =*[“// pu u * d A + u / / pu • dA 

-// P • dA + // (p7u + uVu"^ • dA)] (5) 

To maintain stability, the properties must be updated in the 
proper sequence. In the current method, the sequence is 

1. update the pressure from the continuity equation; 

2. update the velocities from the momentum equation using 
the new pressure and old velocities and old density; 

3. update the density from the ideal gas equation of state; 

4. update the temperature from constant total temperature. 

4. CONTROL VOLUMES 

A new control volume has been introduced for this 
method. To eliminate the need for smoothing of flow proper- 
ties, there must be as many control volumes across the duct 
as there are nodes where these variables are calculated. We 
need as many equations as unknowns. The control volumes also 
need to be located so that errors in continuity and momentum 
can correctly influence the changes in pressure or density 
and velocity without smoothing. The current control volume 
accomplishes this and is shown in Fig. 1. When calculating 
the flux through a streamwise face of an element, the value 
of the flow properties at the node on that face are used. 
Linear Interpolation is used to obtain the flux on the cross- 
stream face. 



Fig. 1 New Control Volumes 
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5. DISTRIBUTION OF PROPERTIES 


The properties at node points are changed In the flow 
field after each time step because the continuity and 
momentum equations are not satisfied for a given control 
volume. A decision must be made about which node, either 
upstream or downstream, these changes should be allocated 
to. The criterion used in determining where changes in 
properties should be sent Is that these distributions result 
In reduced errors In continuity and momentum. The current 
method uses the following allocation procedure: 

1. The pressure is updated through the continuity equation 
and the pressure change is sent to the upstream node; 

2. The u and v velocities are updated through the momentum 
equations and the changes are sent to the downstream 
node; 

3. The density Is updated through the Ideal gas equation of 
state using an Interpolated pressure. 

6. PRESSURE INTERPOLATION PROCEDURE 

As part of the updating procedure used by Denton [5], an 
effective pressure Is used In the momentum equations rather 
than the true thermodynamic pressure determined from the 
Ideal gas equation of state. This effective pressure is 
needed because If the true pressure Is used In the momentum 
equations the solution may not converge. In the current 
method, the density used In the continuity and momentum 
equations Is an effective density which may be different than 
the density obtained using the Ideal gas equation of state. 
This effective density Is used to satisfy stability 
requirements. 

Starting with a generalized pressure Interpolation 
equation for the effective density 


' Pt-1> 

Pl+1 * t^l ®0^^l+l " ^l^ 2 

^^I+l " ^1-2^1 1 

«2 3 ^ » 

Mach number limitations were sought for Sq, aj, and ^2 
that 


( 6 ) 

such 


3i + 32 + 83 “ 1 (7) 

which assures second order accurate solutions. A set of 
equations for sq, a^, and a 2 was chosen which satisfies two 
stability criteria [6]. The equations are 
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for M < 2 


; ai = 


^0 “ 



(\-l) 


1 - 3q ; 32 = 0 : 


for M > 2 


0 ; a. 




1 - ai 


( 8 ) 


These Mach number dependent formulations for aQ, a^, and 82 
are shown In Fig. 2. 



MACH MUMBER 

Fig. 2 Mach Number Dependent Values for 
Coefficients 3 q, ap and a 2 


7. TIME STEPS 

A unique feature of this method Is the use of different 
time steps for the continuity and momentum equations. 
Previous workers who have used explicit time marching methods 
have used the CFL condition as a basis for determining 
allowable time steps which maintain stability. The same time 
step Is used for both the continuity and momentum 
equations. In the current method, the expressions that are 
used to determine the allowable time steps are; for the 
momentum equations 


u 


’'eff 

+ 1 1 

Tk 


5y 

U(6y)^' 


(9) 


and for continuity, 


5t < 
c 



1 



^ 1 a- 

V 


my 


( 10 ) 


where <St is the momentum time step, <S t^ Is the continuity 
time step and is an effective y-coraponent of velocity. 

The advantage of using different time steps is that, for low 
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velocity regions of the £low» the allowable momentum time 
step can be significantly larger than that allowed by the CFL 
condition. These larger time steps allow the boundary layer 
profiles to change more rapidly and enhance the convergence 
rate significantly compared with a method which uses the CFL 
condition. 

8. BOUNDARY CONDITIONS 


For viscous flow, at the upstream boundary, the total 
temperature, freestream total pressure, inlet boundary layer 
velocity profile, and flow angle are specified. Along the 
downstream, boundary the static pressure is specified. Pres- 
sures along the solid boundaries are determined from linear 
extrapolation. For viscous flow, the values of the x- 
component and y-component of velocity are set equal to zero 
at solid walls. 


9. TURBULENCE MODEL 


A Prandtl mixing length model is used to model the 
turbulent stresses. The model is 


‘^eff 


+ M, 



"du" 


L is smaller of 0.08 times the width of boundary layer 
or 0.41 times the distance to the wall 


Van Driest Correction 

L » 0.41 "y"(l - expl- "y” ^/26 Uj^]) 
Near Wall Correction *^eff * ^t^ 


10. MULTI-VOLUME METHOD FOR PRESSURE CHANGES 

Control volumes are grouped in the boundary layer to 

form a larger global control volume. The continuity error is 

calculated for this global control volume and changes in 

pressure are assigned equally to each of the upstream nodes 
for each control volume making up the global control 
volume. Then the global control volume is made successively 
smaller towards the wall. This is shown schematically in 

Fig. 3. The entire pressure change for one iteration at each 
node within the multi-volume region is determined by adding 
together all the pressure changes assigned to that node. 

The multi-volume method propagates pressure changes 
rapidly through the boundary layer and minimizes transverse 
pressure gradients in the intermediate solution. The above 
changes allow the calculation of boundary layer flows where 
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the control volumes near the wall can have aspect ratios 
(length/height) over 1000. 



VOL. I VOL. 2 '^OL. 3 Vql. 4 

Fig. 3 Multi-Volume Method for Pressure Changes 
in the Boundary Layer 

n. TRANSVERSE UPWIND DIFFERENCING 

When the control volumes become long and thin near the 
wall of the duct, the fluxes through the top and bottom faces 
of the control volume become more significant in comparison 
to the fluxes through the streamwise faces. To strengthen 
the diagonal dominance of the coefficient matrix, the 
momentum fluxes through the transverse faces may be calcu- 
lated using interpolated velocities upstream in the 
transverse direction rather than the actual interpolated 
values. The interpolation functions and the derivation of 
the functions is discussed in more detail in Ref. 6. 

12. SAMUEL AND JOUBERT INCOMPRESSIBLE TURBULENT BOUNDARY 
LAYER 

Incompressible turbulent boundary layer flow in a 
diverging duct was calculated for test case 0141 of the 

Stanford Conference [7]. The grid used in the present 
calculations is shown in Fig. 4. The inlet velocity is 26 
m/s. 



Fig. 4 Geometry and Grid for Samuel and Joubert 

Figure 5a shows a comparison of the calculated skin 
friction coefficient with the experimental results and with 
the results from the Moore cascade flow program. The 
agreement is excellent. A comparison of the calculated 
turbulent shear stress distribution, uv, with the experl- 
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mental results is shown in Fig. 5b. The agreement is good. 
Figure 5c shows good agreement also between the calculated 
and measured velocity profiles at two locations in the duct. 

0.000 
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0 002 
0.000 

5a) Skin Friction Coefficient 
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0 00 
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5c) Velocity Profiles 
Fig. 5 Results for Samuel and Joubert 


PLOT 4 CASE 0141 FILES 14.16 



0 0.5 1 0 05 1 

U/U. 



PLOT 1 CASE 0141 FILE 4 



0 0.5 I 1.5 2 2.5 


X (m) 


- 66 - 






13. MDRL DIFFUSER CALCULATIONS 


The diffuser geometry (Model G) Is shown In Figure 6a 
[8,9]. Figure 6a also shows the computational grid used 
which has 87 grid points in the axial direction and 20 points 
across the flow. The Inlet boundary layer thicknesses were 
specified as 9% and 4.5% of the inlet diffuser height for the 
curved and flat wall boundary layers, respectively. For this 
calculation, the ratio of exit static pressure to the inlet 
total pressure was 0.826. In the experiment, this test point 
results in transonic flow In the diverging portion of the 
duct with a Mach number of approximately 1.235 upstream of a 
nearly normal shock, and the flow remained fully-attached 
throughout the diffuser at this test condition. 

A contour plot of static pressure Is shown In Fig. 6b. 
The shock can be seen In the diverging portion of the duct. 
The shock Is well defined as Illustrated by the high 
clustering of contours at the shock. Figure 6c shows a Mach 
number contour plot for the calculations. The calculated and 
measured curved wall static pressures are compared In Fig. 
7. The shock Is well defined and no overshoot occurs In the 



a] 

geometry and 

grid 




b) Static pressure contours 


pMTTlJ 

11* 


c) Mach number contours 


Fig. 6 Geometry and Contours for MDRL Diffuser 

Measured shock locations on the curved wall and In the 
middle of the duct are plotted In Fig. 8 as a function of 
shock Mach number, determined from the minimum wall 

static pressures on the curved wall. The minimum wall static 
pressure In the calculation Is located at x/h » 1.5; this Is 
taken to be the location of the shock. The Mach number 
upstream of the shock was determined to be 1.256 from the 
calculated total pressure ratio across the shock In the 
frees tream. This result Is plotted In Fig. 8 and It agrees 
well with the measured shock location. Comparisons of 
calculated and measured velocity profiles (see Ref. 9) at two 
axial locations along the duct are shown In Fig. 9. The 
agreement Is good. The mass averaged total pressure at the 
diffuser exit divided by the Inlet freestream total pressure 
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A EXPERIMENT 
♦ CALCULATION 
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Fig. 7 Curved Wall Static Pressures for MDRL Diffuser 



CALCULATED 


Fig. 8 Comparison of Computed and Measured 
Shock Position in MDRL Diffuser 



CALCULATION 
A EXPERIMENT 
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+ CALCULATION 
A EXPERIMENT 


Fig. 9 Velocity Profiles at x/h= 4.03 and 8.2 
in MDRL Diffuser 



is calculated from the numerical results to be 0.9615. This 
compares well with the measured value of 0.965, obtained from 
the data of M. Sajben and T. J. Bogar, midway between the 
diffuser side walls. 

The total CPU time for the MDRL diffuser calculations 
was approximately 35 minutes on an IBM 3031. 

14 . CONCLUSIONS 

An explicit finite volume time marching method has been 
extended to allow the calculation of laminar and turbulent 
flow in ducts. Both subsonic and supersonic flow can be 
calculated with the method. Incompressible turbulent 
boundary layer flow in an adverse pressure gradient was 
calculated. The agreement between the calculated and 
measured skin friction coefficient, turbulent shear stress 
distribution, and mean velocity profiles was good. Transonic 
viscous flow through a converging diverging nozzle was 
calculated. The computed and measured velocity profiles were 
In good agreement especially near the exit of the nozzle. 
The computed and measured shock locations were compared and 
were found to be in good agreement. Viscous and shock losses 
in the diffuser were well modelled. 
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APPENDIX C 
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Explicit Finite-Volume Tlme-Marchlng Calculations 
of Total Temperature Distributions In Turbulent Flow 

Stephen Nicholson, Joan G. Moore and John Moore 

Mechanical Engineering Department 
Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 24061 

1 . SUMMARY 

A method has been developed which calculates two-dimen- 
sional, transonic, viscous flow In ducts. The finite volume, 
time marching formulation Is used to obtain steady flow solu- 
tions of the Reynolds-averaged form of the Navler Stokes 
equations. The entire calculation Is performed In the physi- 
cal domain. This paper Investigates the Introduction of a 
new formulation of the energy equation which gives Improved 
transient behavior as the calculation converges. The effect 
of variable Prandtl number on the total temperature distribu- 
tion through the boundary layer Is also Investigated. 

A turbulent boundary layer In an adverse pressure gradi- 
ent (M =« 0.55) Is used to demonstrate the Improved transient 
temperature distribution obtained when the new formulation of 
the energy equation Is used. A flat plate turbulent boundary 
layer with a supersonic frees tream Mach number of 2.8 Is used 
to Investigate the effect of Prandtl number on the dis- 
tribution of properties through the boundary layer. The 
computed total temperature distribution and recovery factor 
agree well with the measurements when a variable Prandtl 
number Is used through the boundary layer. 
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2. INTRODUCTION 


This paper is an extension of the work reported else- 
where In this conference [1], A review of the features of 
the new method will be Included here but a more complete 
discussion may be found In references 1 and 2. 

The features of the current method can be summarized as 
follows. Control volumes are chosen so that smoothing of 
flow properties, typically required for stability. Is not 
needed. Different time steps are used in the different gov- 
erning equations to Improve the convergence speed of the 
viscous calculations. A multi-volume method for pressure 
changes In the boundary layer allows calculations which use 
very long and thin control volumes (length/height = 1000). 

3. GOVERNING EQUATIONS 

The unsteady forms of the continuity equation, the x- 
momentum equation, the y-momentum equation, and the energy 
equation. In Integral form, are used to obtain steady-state 
solutions for flow through 2-dlmensional ducts. This ap- 
proach differs from our previous work [1] where the assump- 
tion of constant total temperature was used Instead of the 
full energy equation. The ideal gas equation of state and a 
Prandtl mixing length turbulence model (1) complete the 
governing equations needed to solve for the unknown vari- 
ables p,u,v,P,M, and T. 

For a finite control volume where we can assign one 
value of density to the control volume, and for a finite time 
step, 6t, continuity states that, 

- p*^ » 5p » -[ JJ p u • d a] (1) 

where the Integral Is evaluated explicitly at the current 
time step, n. In arriving at an expression which relates the 
pressure change directly to the continuity error, we will 
assume that changes in temperature are small In comparison to 
other changes for one time step. Thus, we can relate changes 
In pressure to changes In density through the Ideal gas equa- 
tion of state. 

pn+1 _ pH , 5p , // p u . d a] (2) 

For the method Introduced In the current work, a non-conserv- 
ative form of the unsteady momentum equation Is used. The 
non-conservative form Is used because It allows the current 
method to use different time steps for the continuity, momen- 
tum, and energy equations. The difference between the non- 
conservative and conservative forms of the unsteady momentum 
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equation is associated with the unsteady and convective 
terms. Specifically, we note that 

3(pjj) 3u 

+ V • Puu = p-j^ + pji* 7 u (3) 

and the right hand side of Eq. 3 can be rewritten as 
du 3u 

p^+pu*7u“p^+V»puji-u(7»pu) (4) 

When the right hand side of Eq. A is combined with the pres- 
sure and viscous terms, the momentum equation in integral 
form becomes 

“ 5(ti) « [”// Puu* d A + u ff p u* dA 


(5) 

- ff • d A + // (u 7 u + p V u^) • d A ] 

To maintain stability, the properties must be updated in the 
proper sequence. In the current method, the sequence is: 

1. update the pressure from the continuity equation; 

2. update the velocities from the momentum equations using 
the new pressure and old velocities and old density; 

3. update the density from the ideal gas equation of state; 

4. update the temperature from the energy equation. 

4. ENERGY EQUATION 

For many calculations of transonic viscous flow, the 
assumption of constant total temperature will give a suffi- 
cient representation of the energy equation in the flow 
field. By assuming constant total temperature, the computa- 
tions are less expensive to run and the computer storage 
requirements are less. The assumption of constant total 
temperature is usually satisfactory if: 

1. an adiabatic wall is assumed in the calculations; 

2. no work is done on the fluid at the solid boundaries; 

3. the Mach numbers in the flow fields are low enough that 
total temperature gradients within the boundary layer are 
small; 

4. the Prandtl number is approximately 1.0. 

For a Prandtl number of 0.9, the solution should not 
deviate greatly from the constant total temperature assump- 
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tlon. However for high speed flow, the energy equation 
should be Included in the calculations especially If the 
Prandtl number deviates greatly from 1. 


Two forms of the Integral formulation of the energy 
equation will be derived next. 


The energy equation In differential form is 

3E = 

-y- + 7 • E^u*-7 »^ + 7 • [u» (ii7jj + jj7 jk] - 7 


where the total energy per unit volume, E^, Is 


P (e + V2 (u^ + v^)) = pe^ 


The left hand, side of Eq. 6 can be rewritten as 
3E 9(pe ) 

y^ + 7 . E^u »-y~ + 7 • P e^ u 


and 


3(pej.) 

“TT" 


+ 7 • (pe^) u^ 


^®t 

TT 


+ p u • 7e 


P u 

( 6 ) 

(7) 

( 8 ) 

(9) 


then, expanding the right hand side of Eq. 9, we get, 

p yy- + p ji • 7e^ = p y^ + 7 U e^ - e^(7 • p ji) (10) 

The procedure just outlined Is identical to what was 
done to the unsteady and convective terms In the momentum 
equation (see Eqs. 3,4). 


The heat flux vector, _q^, can be represented as 

q = -k7T 


Substituting Eqs. 8-11 Into Eq, 6, we get 

p yy- * -7 • p u e^ + e^(7 • p ji) — 7 


+ 7 


[2i*(li7u + u7u )]-7 


(-k7T) 
P u 


( 11 ) 


( 12 ) 


The integral form of the energy equation Is then 

6e. 


t 

TT 


p X 6Vol = 


- // p u e^ 


A + f I P u • d A - J f -k7T • d A 


(13) 
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+ // [u« (yVu + yVu'^) »dA-// Pu* dA 

where e^, is an average value for the control volume. As 
with ^ the momentum equation, Eq. 13 has a 
term e^^ // p • d A , which removes the continuity error 
contriDution ”to the^energy error. 

This foirm of the energy equation, when incorporated into 
the current method, behaved poorly. Initially there were 
large errors in continuity and momentum and these large er- 
rors acted through this energy equation to cause errors in 
the total energy for a control volume. This interaction was 
destabilizing. 

An alternative form of the energy equation will now be 
derived. This alternative form has enhanced convergence 
properties when compared with the above formulation. Brief- 
ly, the energy equation is reformulated so that changes in 
total enthalpy, h^., are calculated rather than changes in 
total energy, e^, which was done previously. This allows us 
to see the terms which cause departures from uniform total 
temperature - for both the steady state solution and the 
transient solution. 


The total enthalpy can be defined in terms of the total 
energy and the static temperature 


or 


+ P/p 


e^ + RT 


(14) 

(15) 


Taking the derivative with respect to time and multiplying by 
the density, we get 

3h. 3e, 

V w |- 

(16) 


t 

"TT 


3T 


P Tt" Tt 


The static temperature T can be represented in terms of the 
total enthalpy and the absolute velocity as 


Therefore 



"t 

(17) 


T » — 

TT 2C 
P P 

3T 

^ 1 ^ £V 

(18) 

3t 

“ C 3t ~ C 3t 
P P 

18 into 

Eq. 16, we obtain 


^®t 

P . R „ 3V 

(19) 

^ TT 

’ - IT ^ TT Tt 

P 
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where y Is the ratio of specific heat capacities and V is the 
magnitude of the velocity vector. Using equations (19) and 
(14) to eliminate e^ from equation (12) we get 
3h 

y Tt" ** ^ 

( 20 ) 

+ V«[u • (u7u + WV u"^] - p V II 

" P 

Using h^ » C T + V^/2 and k * u C /Pr, kVT may be replaced 
by P 

2 

kVT - ^ Vh^ - V(|-) (21) 

and from continuity we may replace V»p u with -3p/3t. 

Therefore the energy equation written as a conservation equa- 
tion for total enthalpy is 
3h 

y Tt" ~ ”^*P - ^^t 

(I) (II) (III) (22) 

+ V-y{l - ^)V(^) + V-uCu-V) „ . - gS V 

P 

(IV) (V) (VI) (VII) 

Terms I and II when combined give - p u • 7h^. Therefore 
terms I + II and III contain h^ only in tRe form Vh . Thus, 
when these are the only Important terms in the equation, flow 
with uniform total temperature at the inlet will retain this 
uniform total temperature provided that the boundairy condi- 
tions are consistent with this. 

Term IV is a viscous transport term for total enthalpy when 
the Prandtl number is other than 1. Term V is another vis- 
cous transport term. It however contains the expression 
(u • 7) u which is the gradient of the velocity in the dlrec- 
tTon of *^e velocity; these gradients are usually small com- 
pared with other velocity gradients. Since terms IV and V 
have the form 7 • ( ), they are not source terms, rather they 
can only redistribute the total enthalpy. Terms VI and VII 
on the other hand have the form of source terms. Relative to 
the steady state, they are proportional to the continuity 
error and the momentum error, respectively. We may write 
them as 


M 


I 


3P 

n 


^ 3V 
” Tt 


At the steady state, Eq, 22 becomes 


(23) 
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0 


(24) 


2 

- -V‘p u h^. + 7 . 7h^ + Vp(l - ) 

+ V • u(u • V) u 

Therefore we may artlbrarlly alter the variables 1 and m in 
Eq, 23 and the steady form of the energy equation, Eq. 24, 
will be obtained for converged solutions. The transient 
behavior of h^ is Improved in the calculation procedure by 
choosing 1 » m • 0, l.e, by omitting the transient source 
terms In the enthalpy equation. 

In integral form then the equation for enthalpy changes is 
fiht 

P 6Vol =» y{- //puh^» dA+h^// pji* dA 


(25) 

+ // ^ 7h^ d A + // [(u - u» 7u^ + UU* * <1A 


where u =» u. + 


and ^ 




TrZ 


The time step used for the enthalpy equation is the same as 
for the momentum equation. If the transient source 
P 3p 

term — ^ had been retained in the enthalpy equation, it 

would have been necessary to link the continuity and energy 
equation time steps. Omitting this term allows us to use 
different time steps for the energy equation. 


5. TEST CASES 


Two test cases will be used to explore various aspects 
of the more complete form of the energy equation, Eq. 25, 
discussed previously. 

5.1 Turbulent Boundary Layer in an Adverse Pressure Gradient 

The geometry and grid used in this test case are shown 
in Fig. 1. Flow in this geometry was used in Ref. 1 to test 
the accuracy of the new computational scheme. In Ref. 1 the 
velocities in the duct were low enough that the flow could be 
treated as incompressible. Here, the inlet freestreara Mach 
number was increased to 0.55. The purpose of this test case 
was to illustrate the advantage of the new formulation of the 
energy equation. 

The static temperatures presented in Fig. 2 are from 
calculations after 500 Iterations. It can be clearly seen 
that the new formulation, Eq. 25, gives a better transient 
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solution to the energy equation and It should result In a 
reduction In the computer time required to reach a steady 
state solution* Fig, 3 shows the corresponding total temper- 
ature profiles for the two formulations of the energy equa- 
tion. 



Fig. 1 Grid and Geometry Used to Demonstrate the 
Advantages of the New Formulation of the 
Energy Equation 



Fig. 2 Static Temperature Distribution Through 
the Boundary Layer at M=0.55, x=200 mm, 
after 500 iterations 

5.2 Flat Plate Turbulent Boundary Layer at M • 2.8 

Van Driest [3] presents the total temperature distribu- 
tion within a flat plate turbulent boundary layer with a 
freestream Mach number of 2.8, The experimental total tem- 
perature distribution Is shown in Fig, 4, The geometry and 
grid for these calculations are shown in Fig, 5, The height 
of the duct was 63,5 mm and the length of the duct was 254 
mm. The computational grid shown In Fig, 5 consists of 21 
axial grid points and 14 transverse grid points. The Inlet 
boundary layer thickness of 6,35 ram was 10% of the duct 


- 78 - 


l.Ol T 



l.O 



O 


^ OLD 


0.99 L 
0.0 


0.1 


0.2 


Y/H 

Fig. 3 Total Temperature Distribution Through the 
Boundary Layer at M=0.55, x=200 mm, 
after 500 iterations 

height. The Reynolds number based upon axial distance Is 
approximately 10', To stabilize these supersonic flow 
calculations, the upwind effective density method was used 
[2]. This means that an effective density used at a grid 
point is calculated with the Ideal gas equation of state 
using the pressure from the next upstream grid point. The 
inlet velocity, total temperature, and total pressure were 
specified at the upstream boundary. Three calculations were 
performed with different assumptions about the turbulent 
Prandtl number. These assumptions were 

1. Pr^ = 0.90 Pr. = 0.73 

2. Pr^ =■ 0.73 PrJ = 0.73 

3. Pr,. varies linearly through the boundary layer from 0.9 
at the wall to 0.66 in the frees tream. 


The turbulent Prandtl number Is typically set equal to a 
constant of 0.9 In calculations [4]. The calculated total 
temperature distribution through this boundary layer using a 
constant turbulent Prandtl number of 0.9 is shown in Fig. 6 
(represented as Q ) . The recovery factor is calcu- 

lated to be 0.920 which compares with the empirically deter- 
mined value of 0.90, However, the distribution of total 
temperature through the boundary layer does not compare well 
with the experiment. If the turbulent Prandtl number is set 
equal to the laminar Prandtl number of 0.73, the total tem- 
perature distribution changes as seen in Fig, 6 (represented 
atO's). The distribution through the outer part of the 
boundary layer has Improved but the recovery factor of 0.813 
does not compare well with the experimental value of 0.90. 
Schlichting [5] notes that the turbulent Prandtl number is 
not constant through the boundary layer. The experiments of 
H, Ludwieg [6] for turbulent flow through a pipe show that 
the Prandtl number varies from approximately 0.9 at the pipe 
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Fig. 4 Experimental Total Temperature Distribution 
in a Flat Plate Turbulent Boundary Layer 
M=2.8, Re = 10^ after van Driest 



Fig. 5 Geometry and Grid For Boundary Layer Calculations 
at M=2.8 



Fig. 6 Total Temperature Distribution For Flat Plate 
Boundary Layer at M=2.8 




wall to 0.66 at the center of the pipe. This distribution is 
shown in Fig. 7. The variation is almost linear. For the 
third set of calculations, the Prandtl number was assumed to 
vary linearly through the boundary layer from 0.9 at the wall 
to 0.66 at the edge of the boundary layer. The total temper- 
ature distribution for this case is shown in Fig. 6 (repre- 
sented asA's). The total temperature distribution calcu- 
lated using a variable Prandtl number is also compared with 
the e 3 q>erlmental results in Fig. 8. Both the distribution of 
total temperature through the boundary layer and the recovery 
factor of 0.90 are in good agreement with the experimentally 
measured values. 



Fig. 7 Ratio of the Turbulent Transfer Coefficient Over 
the Length of a Radius in Turbulent Pipe Flow 



Fig. 8 Total Temperature Distribution For Flat Plate 

Boundary Layer at M=2.8 Computation vs. Experiment 
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6. CONCLUSIONS 


A new formulation for the energy equation was Introduced 
which has improved transient behavior when compared with the 
standard formulation. The new formulation removes the in- 
fluences of continuity and momentum errors from the energy 
equation during transients in the solution. 

For flat plate turbulent boundary layer flow with a 
freestream Mach number of 2.8, the calculated total tempera- 
ture profile was Improved by using a variable Prandtl number 
through boundary layer. The recovery factor of 0.90 agreed 
very well with the empirically determined value of 0.9. 
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